A novel approach using nonlinear surfaces for dynamic aperture optimization in MBA synchrotron light sources

MBA cell-based synchrotron light sources have enabled an unprecedented increase in beam coherence and brightness, greatly benefiting the scientific disciplines that rely on X-ray techniques. However, controlling the electron dynamics is a theoretical and technological challenge, due to the large number of parameters to adjust and constraints to satisfy when designing modern synchrotrons. Having versatile tools for the description and manipulation of electron dynamics could favor the design of these accelerators and lead to progress on several fronts in the understanding of matter. In this paper, a formalism based on the use of nonlinear geometric surfaces represented by polynomial quasi-invariants, to analyze and optimize the dynamic aperture of electrons in MBA storage rings, is introduced. The formalism considers on- and off-momentum particle dynamics. Within the optimization scheme, different objective functions defined in terms of the nonlinear surfaces, which are minimized using genetic algorithm methods, are proposed. A remarkable horizontal dynamic aperture exceeding 19 mm is obtained for the design particle of a synchrotron model with 86 pm \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cdot $$\end{document}· rad emittance along with a dynamic aperture above 5 mm for momentum deviations of ± 3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\%$$\end{document}%. According to the results presented, this formalism could be greatly useful for manipulating the dynamical properties of electrons in synchrotrons light sources close to the diffraction limit.

particularly, the use of intense quadrupoles increases the natural chromaticities of the synchrotron 16 .Intense magnetic sextupoles must also be incorporated to bring the chromaticities to zero, with the adverse effect of introducing nonlinearities in the electron dynamics.To reduce these effects, additional sextupoles and higher order magnetic multipoles are incorporated.The combination of multiple constraints and a high number of parameters to be adjusted results in higher demands on the design of the cell's magnetic structure and, of course, technological challenges at various levels.
There are common methods for the optimization of the dynamic aperture of these synchrotrons; for example, the minimization of the resonant terms [17][18][19][20][21] or the direct obtaining of the dynamic aperture by particle tracking [22][23][24][25][26] .However, since the reduction of nonlinear effects in these systems is usually very demanding, several strategies have been proposed and developed, generally based on the aforementioned methods.For example, the combined use of particle tracking and the calculation of resonant terms in optimization objectives 27 ; the proposal of new optimization algorithms or the improvement of existing schemes [28][29][30][31][32][33][34][35][36] ; and the analysis and minimization of fluctuations of conserved quantities in these systems have been explored 28,37,38 .Given the complexity of the problem, these methods usually employ multi-objective functions in which a considerable amount of time is spent to compute their arguments.
In this context, the use of approximate invariant polynomial roots has recently been proposed for the optimization of the nonlinear dynamics of these systems 39,40 .In particular, reference 40 details a dynamic aperture optimization formalism for the design particle, using an objective function dependent on the real roots of the quasi-invariant polynomial of degree 5, since, according to that work, these approximately represent the trajectories of the system.Furthermore, in the above-mentioned study, the dynamic aperture of a storage ring with emittance above 1 nm • rad was optimized showing the advantages of using this type of polynomial in these processes.Whether these polynomials are useful for optimizing electron dynamics in MBA cell-based synchrotrons has not been studied and remained as a pending task.
In this work we present a novel formalism, based on the use of geometric surfaces constructed by quasiinvariants of motion, to study and optimize the dynamic aperture of electrons in a storage ring based on an MBA cell, for both the design particle and considering momentum deviations.This is achieved by extending the studies presented in references [39][40][41][42][43] to propose an approximate polynomial invariant of degree 8 that is used in the construction of nonlinear surfaces, whose contour lines approximately represent the trajectories of the system for a given amplitude and momentum deviation.The concept of nonlinear surface is a fundamental element of the formalism, showing that the construction of quasi-invariants of motion in these systems is a versatile and powerful tool for the analysis and manipulation of their dynamical properties.This versatility is appreciated when, with relative easiness, nonlinear surfaces can be used to define several possible single-objective functions to optimize nonlinear dynamics, allowing a wide search of solutions in the parameter space, which is an advantage over other methods.Within the introduced optimization scheme, these surfaces are incorporated into the proposed singleobjective functions that are optimized using genetic algorithm methods.The use of single-objective functions simplifies the optimization process, and it can be easily done with other numerical optimization methods, which is a significant advantage over other conventional nonlinear optimization schemes 44 .
For the cell design considered, in a configuration with an emittance of 86 pm • rad, a horizontal dynamic aperture greater than 19 mm is obtained for the design particle and greater than 5 mm for momentum deviations of − 3% and 3% .The optimization was performed using a Ryzen-AMD desktop computer in a few hours of computational time.The reduced computer resources needed to perform these calculations is also an advantage over other methods.
These results have been corroborated using the OPA particle tracking module 45 .The optimal nonlinear multipoles were obtained only from the horizontal dynamics optimization.However, when these nonlinear multipoles are introduced in OPA, they generate larger dynamic apertures than when they are optimized with OPA.Based on the results obtained, this formalism could be extremely useful in the manipulation of the dynamical properties of electrons in synchrotrons in the diffraction limit and, therefore, push the frontier of synchrotron development.

Revisiting the Courant-Snyder invariant and approximate constants of motion
In a synchrotron there are several periodic functions associated with the intensities of the different magnetic multipoles, such as the radius of curvature associated with the magnetic dipoles, and the quadrupole, sextupole, and octupole magnetic fields.They are usually represented by ρ(s) , K(s), S(s) and O(s), and satisfy the periodic conditions where c indicates the period, which can be the circumference of the ring, the length of a cell, or a smaller distance inside the cell.As usual 16 , the functions (1) are written in terms of the functions b 1 (s) , b 2 (s) , b 3 (s) and b 4 (s) , defined through the expression (1) , with m = 1, 2, . . ., and Bρ is the magnetic rigidity, which establishes the well-known relationship between the magnetic field, the radius of curvature and the electron energy . Furthermore, for a particle, the deviation of the magnitude of the momentum p with respect to the magnitude of the design momentum p 0 is written as p ≡ p − p 0 and it is proportional to the magnitude of the design moment, �p = δ p 0 , where δ is the percentage of moment deviation with respect to the design moment, i.e. a dimensionless parameter.

Nonlinear surfaces with chromatic effects
Reformulating the discussion found in papers [39][40][41][42][43] with a geometric point of view, the Courant-Snyder invariant is a conserved quantity in the one-dimensional linear system described by the Hamiltonian which allows the definition of the surface which we will call Courant-Snyder surface and describes an elliptic paraboloid whose principal axes depend parametrically on s.Once these are determined for a fixed point s 0 in the ring, the ini- tial conditions (x 0 , p x 0 ) define a particular level curve, associated to a fixed value of the function ( 5): x (which is equivalent to a fixed amplitude of movement).It is well known that the dependence of the principal axes of this elliptic paraboloid on the s parameter (position in the ring) is limited by the periodic conditions imposed on the coefficients appearing in the quadratic function (5) (Courant-Snyder parameters) From this perspective, although the principal axes of the elliptic paraboloid oscillate, the dynamics of the system (4) guarantees that the contour line defined by the equation S 0 (x, p x ) = I 0 , for a given value of I 0 , encloses a constant area as s changes along the trajectory of the ring, despite changing its shape.
In the same way, following the description given in works 39,40,42,43 , the two-dimensional, nonlinear system with chromatic effects, described by the Hamiltonian 17 has the approximate invariant where the functions A (n) ijkl (s) , that generalize the coefficients of the polynomial (3), i.e., α x , β x , γ x , are also subject to periodic conditions Similar to the Courant-Snyder case, we can geometrically interpret the approximate invariant (8) as the hypersurfaces which are parametrically dependent on δ .It is worth noting that, for a fixed value of δ , the surface resulting from the projection of the hypersurface (10) onto the horizontal phase space ( x, p x ) for a fixed value of I has associated contour lines that will not necessarily be simple closed curves.The nonlinearity of the system, introduced by means of the b m (s) coefficients with m > 2 , implies that there will be non-zero f (x, p x , y, p y , δ) , which, as is well known, will result in the aforementioned hypersurface (10) projections (at fixed δ ) having cubic and higher order contour lines, which in general represent unbounded motion.However, studies 39,40 have shown that in the presence of nonlinearities, it is possible to obtain closed contour lines for hypersurface (10) projections (fixed δ ) in horizontal phase space, under certain conditions and for some values of I; i.e., for certain oscillation amplitudes.These curves represent bounded motions of the system and, analogous to the Courant-Snyder case, although the coefficients A (n) ijkl oscillate along the ring, the dynamics of the system (7) allows the area enclosed by each contour line (defined by a fixed value of I) to be approximately conserved.
Due to the context of the problem, we will refer to the projections of the hypersurface (10) onto one of the phase spaces and a fixed δ value, as nonlinear surfaces S δ (x, p x ) .The remainder of this paper will only deal with nonlinear surfaces in horizontal phase space.The fact that the oscillation amplitudes in the horizontal plane are larger than in the vertical one suggests that the dynamics of the system will be more complex in the horizontal than in the vertical plane and, from a physical point of view, it is worth focusing on the horizontal dynamics to advance the understanding of the problem.

Quasi-invariant for a nonlinear, chromatic, two-dimensional system
Following the formalism proposed in works [39][40][41][42][43] , the equations of motion satisfied by the functions A (n) ijkl (s) are obtained by imposing the invariance condition to the expression (8), under the dynamics of the system (7).In this way it is possible to obtain a system of linear differential equations (LDE) whose number of participating functions depends on the degree of the polynomial considered.In order to be able to describe accurately complex structures of resonances in phase space a polynomial of degree 8 is considered in this work, instead of the degree 5 polynomial used in ref. 40 .
A system of 418 equations was obtained using wxMaxima 46 to manipulate this polynomial of degree 8 and at first order in δ .The complete system of equations is not included in this manuscript for ease of reading, but can be consulted in the supplementary material (SM).It should be noted that in this LDE system the equations appear, reproducing the α x , β x and γ x equations, which are important for the consistency of this formalism with the linear case.In addition, the lowest order nonlinear coefficients A   www.nature.com/scientificreports/Only coefficients related to the horizontal dynamics of the accelerator appear in Eqs.(13-14), while in Eq. ( 15) couplings between coefficients associated with both horizontal and vertical space are observed.Within the LDE system, there are many such couplings.On the other hand, the lowest order coefficients of chromatic character A (1) ijkl satisfy the expressions i.e., the coefficients incorporating corrections for the description of off-momentum particles.Note that although these equations contain nonlinear parameters ( b 3 ), in the linear limit ( b n → 0 for n > 2 ) they admit nontrivial solutions.
It is important to consider that these functions will take real and bounded values only if an appropriate selection of b n functions is made.In general, an arbitrary selection of b n functions could produce unstable solutions for the A (n) ijkl functions.Due to the complexity of the 418-dimensional system, finding an analytical solution will require a considerable amount of computational time; therefore, it has been numerically solved in this work imposing the periodic conditions (9).
It is worth outlining the procedure to determine the functions A (0) (s) and the chromatic functions A (1) (s) used in this work.Let us consider a column vector W(s) with 249 components A (0) ijkl (s) followed by 169 components A (1) ijkl (s) , which are determined by the 418 differential equations depicted in SM.A small set of these equations is shown in Eqs.(12-16).
These differential equations can be written in the form where M is a square matrix that depends on dipoles, quadrupoles, sextupoles, and other magnetic multipoles strengths, and can be obtained for each magnetic element k, where the magnetic strengths are considered to be constant.In this case, the general solution of Eq. ( 17), for constant M k , is

) dA
(1) 0011 (1) 0002 Matrices T k , corresponding to the various elements present in the synchrotron cell, can be multiplied in the appropriate order to obtain the transport matrix for the entire cell of p elements, ijkl (s + c) , it can be written representing an eigenvalue problem that can be solved for W(0) if Then, Eq. ( 19) allows to find the values of the periodic functions ijkl (s = 0) contained in the vector W, at s = 0.

Geometric approximation in obtaining bounded motion in phase space
In Ref. 40 it has been shown that it is possible to increase the stable area in the horizontal phase space of electrons in a synchrotron using polynomial quasi-invariants of motion.The proper classification of the real roots obtained from the quasi-invariant polynomial of order 5 into branches poses a serious problem in that work.When this was possible with small oscillation amplitudes, a sextupole optimization could be performed by forcing the chosen branches of the quasi-invariant real roots to topologically approximate the ellipse described by the linear motion through the Courant-Snyder invariant.The procedure for obtaining stable trajectories and resonances is described in Ref. 39 .
In this paper we present a different approach to the above process, i.e., inducing a similarity between the topology of nonlinear and linear phase spaces.To solve the classification of roots difficulty we have opted for the introduction of geometric surfaces described by the quasi-invariant, which seems to be a promising idea according to the results shown below.For this purpose, we will consider the Courant-Snyder (linear) surface S 0 represented by Eq. ( 5).By constructing a grid over the horizontal phase space and calculating the value of the invariant as a function of x and p x on the grid, we obtain this surface as shown in Fig. 1.
The shape of the nonlinear surface, which results from the invariant of Eq. ( 10), depends on the values of the functions A (n) ijkl that contribute the most to the expression of the invariant, and are very sensitive to the configuration of the nonlinear magnetic elements of the cell.Figure 2 is an example of this type of surface, restricted to the horizontal phase space, i.e., y = p y = 0.

Contour lines of nonlinear surfaces
Once the coefficients A (n) ijkl of the polynomial function (10) have been determined for a given value of s, these coefficients can be used to represent nonlinear surfaces for a fixed value of δ .The procedure to achieve this representation will depend on the computational resources available.In this work we have used the surf function of the MATLAB ® numerical computation environment to generate the 3-D surface data and the contour function of the same environment to generate the 2-D contour lines.The contour lines of one of these surfaces, for a fixed value of I, approximately represent the horizontal phase space of the nonlinear chromatic system (7)  to order δ in this work.
As already mentioned, in general, contour lines of nonlinear surfaces do not represent bounded motions in the horizontal phase space, however, in this and the following section protocols to promote the existence of closed contour lines on nonlinear surfaces for different values of δ , under certain conditions, will be described.These protocols are based on optimization schemes that employ the Courant-Snyder surface S 0 (x, p x ) to impose constrictions on the nonlinear surface S δ (x, p x ) , without resorting to the identification of real root branches described in reference 40 , thus avoiding such difficulties.
In this and the next section we consider a grid of dimension n × n , defined over the horizontal phase space, at whose vertices the values of the Courant-Snyder and nonlinear surfaces are computed.Optimizations are studied for moment deviation percentages δ in the range −δ m ≤ δ ≤ δ m , with δ m > 0 the maximum percentage considered.

The proposal of a geometric objective function
In this section two objective function schemes of geometric nature are explored, more robust and simpler to use than the one proposed in Ref. 40 .

Using separation between surfaces
A simple way to define an objective function is to quantify the separation between the surfaces in Figs. 1 and 2 in each grid point and then sum up all these values.This is done for each of the N phase spaces corresponding to the several momenta values δ considered.
This can be done by means of the following expression where S δ (x i , p xj ) and S 0 (x i , p xj ) represents the non linear and linear surfaces evaluated at points (x i , p xj ) of the grid.
(21) For the δ = 0 case the linear surfaces correspond to the Courant-Snyder surface S 0 (x, p x ) , given by Eq. ( 5).To compare a nonlinear surface for δ = 0 with the linear surface for δ = 0 , all the x values of the grid (x, p x ) are displaced to x + D(s)δ , where D(s) is the dispersion function.The nonlinear surface is evaluated using this displaced grid, and then centered at the origin, where the comparison is made.
The definition of f δ obj1 in Eq. ( 21) includes the case for the on-momentum reference particle, i.e, at δ = 0 , as well as terms corresponding to off-momentum particles ( δ = 0 ).These cases will be analyzed separately below.
In minimizing f δ obj1 , we seek for solutions of the nonlinear problem that are as close as possible to the linear surface.As mentioned before, the nonlinear surface is intended to resemble that of the linear case for small amplitudes.
The objective function f δ obj1 defined in the expression (21) considers the sum of distances between the linear and nonlinear surfaces, for each grid point, for a fixed δ phase space.Since these distances depend on the non- linear functions A (0) ijkl (s) and A (1) ijkl (s) , the terms that deform the phase space the most will be penalized the most.Being the only function that is intended to be optimized, the inclusion of weights in the problem is not necessary.We consider this to be a strength of the algorithm by using the objective function of Eq. (21).
Using Gauss curvature A geometrical parameter that can be used as an alternative in the process of optimizing the dynamic aperture in the synchrotron light source is the Gaussian curvature K G .First, it is necessary to know this for a linear surface such as the one shown in Fig. 1.For this purpose, it can be observed that the elliptic paraboloid with semi-axis a, b given by the quadratic form 47 can be approximately expressed as in a neighborhood of the origin (large enough for an MBA lattice optimization), where κ 1 and κ 2 are the principal curvatures.The Gaussian curvature K G is connected with them through Considering the approximation of Eq. ( 23) and that, in a symmetry point of the magnets cell, comparing the expressions (23), (25) and (24), it follows that This last parameter can be a useful guide in the search for a good dynamic aperture in a light source.It can also be useful as a secondary tool after using the objective function f δ obj1 .Linked to this geometrical parameter, a second objective function is defined as where K δ G (x i , p xj ) represents the Gauss curvature surface of the nonlinear surface (Eq. 10),evaluated at points (x i , p xj ) of the considered grid.The curvature K Gauss has been evaluated using the Matlab function surfature 48 through expression where E, F, G are coefficients of the first fundamental form, and L, M, N are coefficients of the second fundamental form.If U ⊂ R 2 , and where h represents a nonlinear surface S δ , then, the coefficients are given by where a subindex letter means derivative with respect to that variable 49 .where β 0 = β x (0) , the Eq. ( 31) can be written in the form where A and B are functions of α , β and γ , provided that θ is chosen in the form to cancel the term x ′ p ′ x .From the expression (34) it also follows that Using (34) and (35) it can be shown that and thus showing that the Gauss curvature of the Courant-Snyder surface is constant, at this level of approximation, at any point s of the synchrotron ring.
In this paper we will refer to the surfaces S δ (x i , p xj ) associated to the quasi-invariant as nonlinear surfaces and in the case of the corresponding Gaussian curvature surfaces K δ G (x i , p xj ) , as associated nonlinear surfaces and, in general, simply as surfaces, when there is no ambiguity.

Application to an MBA synchrotron of the nonlinear-surfaces based objective functions
The use of the methodology described above to achieve good dynamic performance in MBA cell models employed in synchrotron sources in the diffraction limit is studied in this section.

An MBA lattice for the Mexican synchrotron light source combining the ESRF and SLS models, a case of study
In search of an appropriate design for the Mexican light source, a combination of novel aspects, such as the use of LBG 14,50 , of the European ESRF-EBS synchrotron and the design contemplated for the upgrade of the Swiss synchrotron, SLS-2, has been considered.The section that includes a dispersion bump is taken from the European ESRF-EBS synchrotron, and the the internal cells that have been optimized with reverse bends are taken from (31) I = γ x (s 1 )x 2 + 2α x (s 1 )xp x + β x (s 1 )p 2 x .
(32) the Swiss synchrotron, SLS-2.Both models, with emittances of the order of 130 pm • rad, are considered fourth- generation storage rings 51 .A cell of the model under study is shown in Figs. 3 and 4.
The linear optical functions for one cell are shown in Fig. 3.The MBA cell is shown at the bottom of the figure, where dipoles, quadrupoles, and sextupoles are depicted in blue, red, and green, respectively.The complete ring has 20 cells.In Fig. 4 the cell is shown in more detail, depicting the magnetic fields of each magnet, following the magnets color code.
Once the dipole and quadrupole magnetic fields have been fixed, the storage ring optical functions and main parameters are determined; the latter are shown in Table 1.Natural chromaticities indicate that the sextupole intensities should not be very large, therefore, when performing chromaticity corrections, undesirable nonlinear phenomena should be minimal.
This synchrotron model based on an MBA cell represents a suitable problem to study with this methodology due to the complexity of its design, its low emittance and the high number of parameters to be optimized.

Figure 5.
Phase spaces for various values of δ .On the top row, for − 3, − 2, and −1 %.In the center for 0 % and in the lower row, for 1, 2, and 3 %.The optimization was only performed for the phase space corresponding to δ = 0 using the objective function f δ obj1 of Eq. ( 21).The emittance of this operating point is 113 pm • rad.
For this first model, with an emittance of 113 pm • rad, the phase space at δ = 0 has been optimized for several amplitudes, up to x 0 = 7 mm.However, stability is achieved a little beyond this amplitude where it is observed that invariant tori exist bounding resonances of high order.
It is interesting to note that, although the nonlinear surfaces related to the phase spaces for δ = 0 did not intervene in this optimization process, the selected nonlinear parameter scheme also results in zones of stability in these spaces, for x 0 ≤ 5.2 mm in the range of δ ∈ (− 3, 3)%.
To compare with resonant techniques, Table 2 shows the output of the OPA Non-Linear Dynamics module with the set of sextupoles and octupoles that produces the orbits of Fig. 5, without further optimization.Several chromatic and geometric integrals appear with reduced values.Some integrated values of sextupole intensities are close to zero and these sextupoles could probably be discarded.Other values suggest sextupole coupling.The maximum intensity value of sextupoles is not high if we consider that the model has an emittance of 113 pm • rad.It is worth noting that the values of resonant integrals and sextupoles and octupoles integrated strengths are within acceptable ranges even though only horizontal space was considered in the optimization.
Using objective function f δ obj1 : optimizing three phase spaces for δ = − 2, 0, 2 % The idea behind the objective function of Eq. ( 21) is to promote similarity between the topologies of the phase spaces corresponding to nonlinear dynamics, with those of linear dynamics.In this section the range of influence that the functions A (1) ijkl have on dynamic stability is explored.These functions represent the first contribution in δ to the polynomial quasi-invariant.The cell and point of operation are the same as in "Using objective function f δ obj1 : optimizing parameters for the phase space corresponding to δ = 0" section, which has an emittance of 113 pm • rad.
In this case, the optimization process begins with the addition of two functions f δ obj1 , corresponding to offmomentum phase spaces δ = − 2, 2 % , to the objective function f δ=0 obj1 .The amplitude used in the optimization is x 0 = 7 mm.Figure 6 shows the nonlinear surfaces, the trajectories in phase space, and their comparison with tracking simulation, through OPA, generated with the sextupole and octupole configuration obtained with the optimization process using genetic algorithms.A trajectory with amplitude close to 4 mm is shown in this Table 2. Output of the OPA Non-Linear Dynamics module with the set of sextupoles and octupoles that produces the orbits of Fig. 5, without additional optimization.The model has an emittance of 113 pm rad.figure.At δ = 0 (third row) very good overlap is observed between both curves.As δ moves away from 0 while less overlap is observed, the general form of the curve is similar, indicating that the theoretical contribution of A (1) ijkl is less efficient in describing the orbits at large momenta, but it is a good measure of non-linearity.Further study is required to understand whether contributions from the functions A (2) ijkl , which are proportional to δ 2 , are necessary to obtain better overlap at δ away from 0. Having a better theoretical representation of the nonlinear effects should improve the optimal selection of nonlinear elements that increase the dynamic aperture, including off-momentum particles.
In Fig. 7 the most complete panorama is shown from the point of view of particle tracking.Using the set of sextupoles and octupoles found in the optimization process, OPA shows several closed orbits at different amplitudes and momenta.Although the optimization was carried out using only the phase spaces corresponding to δ = − 2, 0 and 2 % in the objective function, higher momenta have been explored with OPA also showing stability, for example, at δ = − 3, and 3 %.The optimization, done with an amplitude x 0 = 7 mm, provides stability up to 4.7 mm at − 3%, 5.9 mm at − 2%, 7.0 mm at − 1%, 8.7 mm at 0%, 10 mm at 1%, 9.7 mm at 2%, 7.5 mm at 3%.By optimizing the 3 phase spaces considered, it is found that there are phase spaces that have slight improvements in the stability zone, compared with one phase space optimization, in exchange for negatively affecting others.These results show that by including the terms A (1) ijkl (s) for δ = 0 , now implemented for treating off-momentum particles, there is stability up to 5.9 mm in the phase spaces for the range of − 2% < δ < 2% here optimized and, it is interesting to note that in the phase spaces of − 3% < δ < 3% there is stability up to 4.7 mm.
The results obtained may be due to multiple factors.Probably the comparison of linear and nonlinear surfaces in the objective function for phase spaces where δ = 0 impose a significant restriction on the topology of the additional spaces, reducing their region of stability.Relaxing the restriction on the similarity of topologies between linear and nonlinear surfaces could lead to more promising results assuming that the A (1) ijkl functions adequately approximate the nonlinear dynamics, behaviour that has been observed in particle tracking simulation for δ ∈ (∼ −2, ∼ 2) .The introduction of contributions of A (2)   ijkl functions could also improve the description of the dynamics in the interval considered and probably extend the range of applicability to values of δ beyond this interval.Research in this direction is currently underway to increase the reliability of the δ-dependent results.

Using objective function f δ
obj2 : optimizing the phase space for δ = 0% In this section the use of the objective function built with the Gaussian curvature, as shown in Eq. ( 27), will be studied.A major problem in synchrotron design is the choice of the appropriate operating point.In the course of this work, a new, lower-emittance operating point emerged, so it was decided to study it in more detail.For this purpose, some linear parameters were changed, including free space lengths.The parameters of this operating point are: the circumference is 484.452m, the emittance is 89 pm rad, the chromaticities are ξ x = − 85.91 and ξ y = − 56.17 ; the tunes are ν x = 48.898, ν y = 26.5386 .Again, we find that the phase advance in the x direction between dispersion bumps satisfy the − I transformation since �ϕ x = 2.9972 π but �ϕ y = 1.9922 π is not a −I transformation 54 .
Figure 8 shows the phase spaces obtained when the optimization is performed with the objective function of Eq. ( 27) that uses the Gauss curvature, at an oscillation amplitude of x 0 = 20 mm, an amplitude that has been Phase spaces for various δ values.On the top row, for − 3, − 2, and − 1 %.In the center for 0 % and in the lower row, for 1, 2, and 3 %.The optimization was performed for the phase spaces corresponding to δ = − 2, 0, and 2 %, using the objective function f δ obj1 of Eq. ( 21).The emittance of this operating point is 113 pm rad.
At this level, we are interested in studying whether the previously obtained dynamic stability can be increased, using this new operating point, by optimizing only the phase space at δ = 0 .For this, the two objective functions f δ obj1 and f δ obj2 will be used.a)Optimization using the f δ obj1 for δ = 0 The solution found during the optimization process of sextupoles and octupoles for this new operating point with 86 pm rad emittance using f δ obj1 for δ = 0 is presented in Fig. 10a.Their maximum stability amplitudes are 4.9, 13.2, − 5.3 mm for δ = −3, 0, 3% respectively.Again, as in Ref. 40 , when the phase space of δ = 0 is optimized, there is stability spillover into phase spaces of δ = 0 .An equivalent stability area is observed for off-momentum phase spaces.There are no significant changes in the stability amplitude for off-momentum phase spaces compared with Fig. 9.   10b.The maximum oscillation ampli- tudes that define the stability regions are 2.1 mm, 21.1 mm and − 5 mm for δ = − 3, 0, 3% respectively.An improvement of the stability region for the on-momentum phase space and a reduction for off-momentum phase spaces can be appreciated when compared with the results obtained with f δ obj1 (Fig. 10a).Both results show a good stability zone for on-momentum phase space and stability spread for the off-momentum phase spaces, for this MBA-based model.This spread has been also observed in Ref. 40 for a third-generation synchrotron model.The off-momentum stability obtained for the MBA-model is rather limited.This probably indicates that the constraint used to achieve linear-like on-momentum dynamics may be too strong.A more relaxed condition may improve off-momentum stability.In the next section, a simple way to relax the condition imposed on the above objective functions is shown.

Relaxing constraints
We consider that it is of interest to ask about the changes that can be produced in the stability of the phase spaces if the restrictions of the objective function are relaxed.In the previous case, the optimization process induces the Gaussian curvature of the nonlinear surface to be close to 4, which is the approximate value of this variable in the linear surface.
A slight decoupling can occur if, in the optimization, each principal curvature is separately required to comply with the relationship in Eq. ( 26), i.e. κ 1 = 2γ x (0) and κ 2 = 2β x (0) in each phase space at fixed δ .Both contribu- tions are added to form the objective function where the superscript 0 indicates that this function is for δ = 0 phase space.
Starting from the solution shown in Fig. 10b, an optimization with this objective function was performed.The phase spaces found are shown in Fig. 11.The maximum amplitudes of the stability zone are 5.8 mm at − 3 %, 19.3 mm at 0 % and − 5.1 mm at 3 %.Good stability zone for each value of momentum δ is obtained and shown in Fig. 12.The zone of greatest stability is close to δ = 0 with amplitudes of horizontal stability x ∼ 20 mm.When the value of δ moves away from 0, the stability zone decreases following a triangle-like geometry.Similar behaviour should be expected in other points around the cell; however, this analysis is beyond the scope of this work.With this scheme there is an increase in the stability zone for the δ = 0 phase spaces compared (41) with the solution shown in Fig. 10b.As a matter of fact, this objective function provides the best solution for the MBA-based model considered.This result shows that reducing the constraints imposed on the objective function allows increasing the stability zone of the off-momentum phase spaces.Schemes of less restrictive objective functions are under study.
Table 3 shows the maximum stability amplitudes of the models presented in the previous sections using the different objective functions in the optimization process.
A more realistic analysis of the quality of the magnetic multipole sets found so far can be achieved with the OPA dynamic aperture calculation module, when using the multipole set optimized with the methods proposed in this work.Figure 13 presents dynamic apertures for some of the models, showing acceptable horizontal amplitudes for y = 0 , despite the low emittances of the models and considering that the vertical dynamics is not optimized.It is worth noting the positive effect on the vertical space, as depicted in Fig. 13c, when the horizontal dynamic optimized set of multipoles employed in Fig. 12 is used.Figures 11 and 13 are calculated with the same set of multipoles optimized with the present theory.Figures 11 and 13 depict 1-D and 2-D results respectively.There is an apparent difference between both figures since the maximum amplitude x ∼ 19 mm of Fig. 11 for δ = 0 is not shown in Fig. 13c for y = 0 .The reason of this apparent discrepancy is that by default DA OPA module performs calculations for y > 0.001 mm and consequently the point (x, y) = (19.3,0) mm does not appear in Fig. 13c.This is also evident in Table 2, where most of the resonant terms related to the vertical space reduce their values.The high values of the tune shift with amplitude dQyy and cubic chromaticity CrY cub terms could be due to the lack of vertical dynamics optimization.Nevertheless, the optimization performed with the method proposed in this paper generates a higher dynamic aperture than that obtained by using the OPA Nolinear Optimization module.
Additionally, the behavior of the maximum stability amplitude as the objective function reduces its value in the optimization process is presented.Figure 14 shows, for a particular optimization case, that the maximum stability amplitude increases as the normalized objective function takes smaller values.The horizontal axis is taken as the smaller of the maximum stability amplitude |x|.The relation between both variables is non-smooth, probably due to the sign ambiguity of x taken at p x = 0 , and the fact that the objective function quantifies the difference between linear and non-linear surfaces.Although the dynamic aperture is not directly optimized, there is a clear relation between the value of the objective function and the maximum stability amplitude.

Concluding remarks
A novel method based on nonlinear surfaces and surface curvatures, derived from quasi-invariant polynomials, has been successfully applied to optimize the nonlinear dynamics of a synchrotron MBA model.The primary objective was to identify maximum amplitude phase space trajectories.The complexity of the chosen MBA model, its low emittance and the high number of parameters to be optimized, make it suitable to assess method's applicability.This evaluation suggest its robustness and suitability to optimize nonlinear dynamics in modern synchrotron light sources.The study has been extended to other lower emittance lattices not reported in this Table 3. Summary of the main parameters and values of the maximum stable amplitudes of the studied models, obtained with the objective functions f 1 = f 0 obj1 , f 2 = f 0 obj2 , f 3 = f 0 obj3 , f δ = δ∈{−2,0,2%} f δ obj1 , for three phase spaces ( δ = −3, 0, 3%).paper, showing equivalent results to those of this work, suggesting that an adequate description of the phase space is provided by the 8th degree polynomial quasi-invariant.
The control of nonlinear phenomena that appears in the studied synchrotron optimization, with a large number of free parameters (20 sextupoles and 6 octupoles), has been considered in order to test the tools used in the work, with encouraging results.
The optimization, based on the construction of a quasi-invariant polynomial of degree eight, considers that the nonlinear surfaces introduced in this work should be topologically similar to those of Courant Snyder or linear surfaces.
It has been demonstrated, at multiple operating points of the synchrotron model, that optimizing the onmomentum phase space leads to a wide stability region that extends to off-momentum phase spaces.The stability zone reduces as the momentum dispersion deviates from zero.Due to the broad on-momentum stability zone achieved, further off-momentum optimizations may benefit using the obtained nonlinear magnetic parameters as an starting point.
The nonlinear surfaces proportional to the momentum dispersion has been introduced, in order to adequately weigh off-momentum spaces.The results suggest that the description using A (1) ijkl functions may be insufficient for |δ| 2% and, probably, the use of A (2) ijkl functions will be needed in order to properly describe this off-momentum range.
The improved results obtained with the introduction of the more relaxed function f 0 obj3 suggest that the constraint on the objective function from similarity between nonlinear and linear surfaces may be too restrictive and this condition could be relaxed.Work in this direction is ongoing.

Figure 1 .
Figure 1.Linear surface corresponding to the Courant-Snyder invariant at the start of the unit cell.Colors correspond to different values of the invariant S 0 (x, p x ) .Axes units are [m] and [rad], corresponding to the coordinate x and conjugate momentum p x of the horizontal phase space, and [m rad] for S 0 (x, p x ).

Figure 2 .
Figure 2. (a) Representative nonlinear surface constructed from a polynomial quasi-invariant of degree 8 corresponding to a phase space for a momentum deviation δ . (b) Contour lines of the previous surface. (c) Phase space of the same system obtained with the OPA particle tracking module.

Figure 3 .
Figure 3. Optical functions of an MBA cell in the synchrotron ring considered.The functions β x (s) and β y (s) are respectively marked in blue and red and the dispersion function η(s) is marked in green.The lower part shows the distribution of dipoles (blue), quadrupoles (red), sextupoles (green); and octupoles are represented with brown vertical lines of the MBA cell.

Figure 4 .
Figure 4. Magnetic fields of the magnets in Fig. 3 are shown.The left set of magnets, separated by drift spaces, shown in the bottom of the figure, corresponds to a matching cell.The field strengths for the corresponding dipoles (blue), quadrupoles (red), and sextupoles (green) are depicted in the upper part of the figure.These field values are related with the parameters b n through Eq. (2).

Figure 6 .
Figure 6.Phase spaces for various δ values.The first column of this figure shows the surfaces S δ (x i , p xj ) corresponding by rows to δ = − 2, − 1, 0, 1, and 2% .The closed paths, confined by the basin, with amplitudes around 4 mm are shown in red.The second column shows the closed paths in phase space and their comparison with the tracking simulation with OPA, in blue.The optimization was performed including the phase spaces corresponding to δ = − 2, 0, and 2 %, using the objective function f δ obj1 of Eq. (21).

Figure 7 .
Figure 7.Phase spaces for various δ values.On the top row, for − 3, − 2, and − 1 %.In the center for 0 % and in the lower row, for 1, 2, and 3 %.The optimization was performed for the phase spaces corresponding to δ = − 2, 0, and 2 %, using the objective function f δ obj1 of Eq. (21).The emittance of this operating point is 113 pm rad.

Figure 9 .
Figure 9. Phase spaces for − 3, − 2, − 1, 0, 1, 2, and 3 % δ values for the new operating point obtained by changing the quadrupoles as indicated in Eq. (39).Sextupoles and octupoles strengths correspond to those used in Fig. 8 without further optimization.The phase spaces for δ = 0 show the stability spread for momenta close to δ = 0 .Different scales are given to figures for clarity.

Figure 10 .
Figure 10.Phase spaces for − 3, 0, and 3 % δ values for the new operating point obtained by changing the quadrupoles as indicated in Eq. (39).Two different solutions were obtained optimizing sextupole and octupole strengths, contrary to what was shown in Fig. 9, with f δ obj1 (a) and with f δ obj2 (b).

Figure 11 .
Figure 11.Phase spaces obtained by optimizing f 0 obj3 for δ = 0 .They are shown in upper row for δ − 3, − 2, − 1 , in middle row for 0, and in bottom row for 1, 2, 3 %.When compared with Fig. 10b, the stability zone for negative and positive momentum deviations has now similar values, increasing the stability amplitude at δ < 0 momenta deviation.

Figure 12 .
Figure12.Space of momentum deviations δ vs horizontal oscillation amplitude x, showing the reduction of the oscillation amplitude as δ moves away from 0. The white area represents the stability zone of the solution shown in Fig.11.

Figure 13 .
Figure 13.Dynamic apertures provided by OPA for δ = 0 for three cases of nonlinear magnetic multipoles optimization, (a) for the case of 113 pm rad of Fig. 5, (b) for the case of 86 pm rad of Fig. 10, (c) for the case of 86 pm rad of Fig. 11.The OPA code does not perform the dynamic aperture evaluation for y < 0.001 mm, thus, data corresponding to y = 0 are not depicted in these figures.

Figure 14 .
Figure 14.This figure shows, for a typical case, the general behavior between the objective function and the maximum stability amplitude.The vertical axis represents the objective function normalized to its initial value, and the horizontal axis represents the smallest |x| in mm, of the maximum stability amplitude, taken at p x = 0.